pacman::p_load(sjPlot, jtools, dplyr, ggplot2, ggpubr, broom, pollster, tidyr, tidyverse, ggthemes, haven)

getwd() # set working directory to folder with replication files in it

setwd("C:/Users/jonsc/Box/ISCAP/Papers/POQ Materials/")

gss <- read_dta("gss7224_r1.dta")

# data cleaning
gss <- gss %>%
  mutate(
    # Party ID simplified (with leaners as partisans)
    pid = case_when(
      partyid %in% c(-99, -98, 7) ~ NA_character_,
      partyid %in% 0:2 ~ "Democrat",
      partyid == 3 ~ "Independent",
      partyid %in% 4:6 ~ "Republican"
    ),
    
    # Party ID with leaners as independents
    partyid_without_leaners_as_partisans = case_when(
      partyid %in% c(-99, -98, 7) ~ NA_character_,
      partyid %in% 0:1 ~ "Democrat",
      partyid == 2 ~ "Independent",
      partyid == 3 ~ "Independent",
      partyid == 4 ~ "Independent",
      partyid %in% 5:6 ~ "Republican"
    ),
    
    # Dummy variables for partisanship
    republican = ifelse(pid == "Republican", 1, 0),
    democrat = ifelse(pid == "Democrat", 1, 0),
    independent = ifelse(pid == "Independent", 1, 0),
    
    republican_noleaners = ifelse(partyid_without_leaners_as_partisans == "Republican", 1, 0),
    democrat_noleaners = ifelse(partyid_without_leaners_as_partisans == "Democrat", 1, 0),
    independent_noleaners = ifelse(partyid_without_leaners_as_partisans == "Independent", 1, 0),
    
    # Age cleanup and category
    age = na_if(age, -100),
    age = na_if(age, -98),
    age = na_if(age, -99),
    age_cat_4 = case_when(
      age < 24.5 ~ "18 to 24",
      age < 44.5 ~ "25 to 44",
      age < 64.5 ~ "45 to 64",
      age >= 64.5 ~ "65 and over"
    ),
    
    # Sex
    sex = recode(as.character(sex),
                 `1` = "Male", `2` = "Female",
                 `-100` = NA_character_, `-99` = NA_character_,
                 `-98` = NA_character_, `-97` = NA_character_),
    sex = factor(sex, levels = c("Female", "Male")),
    
    # Race
    race = recode(as.character(race),
                  `1` = "White", `2` = "Black", `3` = "Other",
                  `-100` = NA_character_),
    
    # Region
    region = recode(as.character(region),
                    `1` = "New England", `2` = "Middle Atlantic",
                    `3` = "East North Central", `4` = "West North Central",
                    `5` = "South Atlantic", `6` = "East South Central",
                    `7` = "Region_West South Central", `8` = "Mountain",
                    `9` = "Pacific"),
    
    # Education & Degree
    education = recode(as.character(degree),
                       `0` = "No College Degree", `1` = "No College Degree",
                       `2` = "No College Degree", `3` = "College Degree", `4` = "College Degree",
                       `-99` = NA_character_, `-98` = NA_character_, `-97` = NA_character_),
    education = factor(education, levels = c("No College Degree", "College Degree")),
    
    degree = recode(as.character(degree),
                    `0` = "Less than High School", `1` = "High School Degree",
                    `2` = "College Degree", `3` = "College Degree", `4` = "Graduate Degree",
                    `-99` = NA_character_, `-98` = NA_character_, `-97` = NA_character_),
    degree = factor(degree, levels = c("Less than High School", "High School Degree", "College Degree", "Graduate Degree")),
    
    # Religion
    religion = recode(as.character(relig),
                      `1` = "Protestant", `2` = "Catholic", `3` = "Other",
                      `4` = "Atheist/Agnostic/No Religion",
                      `5` = "Other", `6` = "Other", `7` = "Other",
                      `8` = "Other", `9` = "Other", `10` = "Other",
                      `11` = "Other", `12` = "Other", `13` = "Other",
                      `-99` = NA_character_, `-80` = NA_character_,
                      `-97` = NA_character_, `-98` = NA_character_),
    religion = factor(religion, levels = c("Protestant", "Other", "Catholic", "Atheist/Agnostic/No Religion")),
    
    # Urban/Rural
    rur_urb = recode(as.character(xnorcsiz),
                     `1` = "Urban", `2` = "Urban",
                     `3` = "Suburban", `4` = "Suburban", `5` = "Suburban", `6` = "Suburban",
                     `7` = "Rural", `8` = "Rural", `9` = "Rural", `10` = "Rural",
                     `-99` = NA_character_, `-100` = NA_character_, `-80` = NA_character_),
    
    rur_urb_binary = case_when(
      rur_urb == "Rural" ~ "Rural",
      rur_urb %in% c("Urban", "Suburban") ~ "Non-Rural"
    ),
    rur_urb_binary = factor(rur_urb_binary, levels = c("Rural", "Non-Rural")),
    
    # Trust in scientists
    consci = recode(consci, `1` = 3, `3` = 1, `-100` = NA_real_,
                    `-98` = NA_real_, `-99` = NA_real_, `-97` = NA_real_),
    consci_binary = ifelse(consci == 3, 1, 0),
    consci_cat = recode(as.character(consci),
                        `1` = "Hardly any", `2` = "Only some", `3` = "A great deal"),
    consci_cat = factor(consci_cat, levels = c("Hardly any", "Only some", "A great deal")),
    consci_cat_binary = ifelse(consci_cat == "A great deal", "Percent 'A great deal' of trust in scientists", 0),
    consci_cat_lowest = ifelse(consci_cat == "Hardly any", "Percent 'hardly any' of trust in scientists", 0),
    
    # Religiosity
    attend = recode(as.character(attend),
                    `8` = "More than once a week", `7` = "Once a week",
                    `6` = "Once a week", `5` = "A few times a month",
                    `4` = "A few times a month", `3` = "A few times a year",
                    `2` = "A few times a year", `1` = "Once a year or less",
                    `0` = "Never", `-100` = NA_character_,
                    `-98` = NA_character_, `-99` = NA_character_, `-97` = NA_character_),
    attend = factor(attend, levels = c("Never", "Once a year or less", "A few times a year",
                                       "A few times a month", "Once a week", "More than once a week")),
    religiosity = case_when(
      attend %in% c("Never", "Once a year or less", "A few times a year", "A few times a month") ~ "Not Religious",
      attend %in% c("Once a week", "More than once a week") ~ "Religious"
    ),
    religiosity = factor(religiosity, levels = c("Religious", "Not Religious")),
    
    # Class
    class3 = case_when(
      class %in% c(1, 2) ~ "Lower/Working",
      class == 3 ~ "Middle",
      class == 4 ~ "Upper"
    ),
    class2 = case_when(
      class %in% c(1, 2) ~ "Lower/Working",
      class == 3 ~ "Middle/Upper",
      class3 == "Upper" ~ "Middle/Upper"
    ),
    
    # Race binary
    race_binary = case_when(
      race == "White" ~ "White",
      race == "Black" ~ "Black"
    )
  )



# Figure 1
crosstab(gss, consci_cat, year, pct_type = "col", weight = wtssps, format = "long") %>%
  mutate(year = as.integer(as.character(year))) %>%
  ggplot(aes(x = year, y = pct / 100, color = consci_cat, group = consci_cat)) +
  geom_smooth() +
  geom_point() +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0, 1)) +
  theme_bw(base_size = 14) +
  xlab("") + 
  ylab("Percentage of Americans") +
  labs(color = "Confidence in\nScientific\nCommunity") +
  scale_color_manual(values = c("gray35", "gray70", "gray0")) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  scale_x_continuous(breaks = c(
    1972, 1973, 1974, 1975, 1976, 1977, 1978, 1980, 1982, 1983, 1984, 1985, 1986, 1987, 1988, 1989, 1990, 1991,
    1993, 1994, 1996, 1998, 2000, 2002, 2004, 2006, 2008, 2010, 2012, 2014, 2016, 2018, 2021, 2022, 2024
  ))


# Table 1
# Calculate mean trust in all institutions by year
gss$conarmy <- car::recode(gss$conarmy, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$conbus <- car::recode(gss$conbus, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$conclerg <- car::recode(gss$conclerg, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$coneduc <- car::recode(gss$coneduc, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$confed <- car::recode(gss$confed, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$confinan <- car::recode(gss$confinan, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$conjudge <- car::recode(gss$conjudge, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$conlabor <- car::recode(gss$conlabor, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$conlegis <- car::recode(gss$conlegis, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$contv <- car::recode(gss$contv, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$confed <- car::recode(gss$confed, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$conmedic <- car::recode(gss$conmedic, '-100=NA; -98=NA; -99=NA; -97=NA')
gss$conpress <- car::recode(gss$conpress, '-100=NA; -98=NA; -99=NA; -97=NA')

gss$conarmy_binary <- ifelse(gss$conarmy==1,1,0)
gss$conbus_binary <- ifelse(gss$conbus==1,1,0)
gss$conclerg_binary <- ifelse(gss$conclerg==1,1,0)
gss$coneduc_binary <- ifelse(gss$coneduc==1,1,0)
gss$confed_binary <- ifelse(gss$confed==1,1,0)
gss$confinan_binary <- ifelse(gss$confinan==1,1,0)
gss$conjudge_binary <- ifelse(gss$conjudge==1,1,0)
gss$conlabor_binary <- ifelse(gss$conlabor==1,1,0)
gss$conlegis_binary <- ifelse(gss$conlegis==1,1,0)
gss$contv_binary <- ifelse(gss$contv==1,1,0)
gss$conmedic_binary <- ifelse(gss$conmedic==1,1,0)
gss$conpress_binary <- ifelse(gss$conpress==1,1,0)

armytrend <- crosstab(gss, conarmy_binary, year, pct_type = "col", weight = wtssps, format = "long")
armytrend %>% filter(conarmy_binary == 1, year == 2024) %>% pull(pct)

scientisttrend <- crosstab(gss, consci_binary, year, pct_type = "col", weight = wtssps, format = "long")
scientisttrend %>% filter(consci_binary == 1, year == 2024) %>% pull(pct)

bustrend <- crosstab(gss, conbus_binary, year, pct_type = "col", weight = wtssps, format = "long")
bustrend %>% filter(conbus_binary == 1, year == 2024) %>% pull(pct)

clergtrend <- crosstab(gss, conclerg_binary, year, pct_type = "col", weight = wtssps, format = "long")
clergtrend %>% filter(conclerg_binary == 1, year == 2024) %>% pull(pct)

eductrend <- crosstab(gss, coneduc_binary, year, pct_type = "col", weight = wtssps, format = "long")
eductrend %>% filter(coneduc_binary == 1, year == 2024) %>% pull(pct)

fedtrend <- crosstab(gss, confed_binary, year, pct_type = "col", weight = wtssps, format = "long")
fedtrend %>% filter(confed_binary == 1, year == 2024) %>% pull(pct)

finantrend <- crosstab(gss, confinan_binary, year, pct_type = "col", weight = wtssps, format = "long")
finantrend %>% filter(confinan_binary == 1, year == 2024) %>% pull(pct)

judgetrend <- crosstab(gss, conjudge_binary, year, pct_type = "col", weight = wtssps, format = "long")
judgetrend %>% filter(conjudge_binary == 1, year == 2024) %>% pull(pct)

labortrend <- crosstab(gss, conlabor_binary, year, pct_type = "col", weight = wtssps, format = "long")
labortrend %>% filter(conlabor_binary == 1, year == 2024) %>% pull(pct)

legistrend <- crosstab(gss, conlegis_binary, year, pct_type = "col", weight = wtssps, format = "long")
legistrend %>% filter(conlegis_binary == 1, year == 2024) %>% pull(pct)

medictrend <- crosstab(gss, conmedic_binary, year, pct_type = "col", weight = wtssps, format = "long")
medictrend %>% filter(conmedic_binary == 1, year == 2024) %>% pull(pct)

presstrend <- crosstab(gss, conpress_binary, year, pct_type = "col", weight = wtssps, format = "long")
presstrend %>% filter(conpress_binary == 1, year == 2024) %>% pull(pct)

tvtrend <- crosstab(gss, contv_binary, year, pct_type = "col", weight = wtssps, format = "long")
tvtrend %>% filter(contv_binary == 1, year == 2024) %>% pull(pct)


# Figure 2
gss$consci_gender <- ifelse(gss$consci_cat_binary=="Percent 'A great deal' of trust in scientists", "Gender", 0)
gss$consci_race <- ifelse(gss$consci_cat_binary=="Percent 'A great deal' of trust in scientists", "Race", 0)
gss$consci_urb <- ifelse(gss$consci_cat_binary=="Percent 'A great deal' of trust in scientists", "Rurality", 0)
gss$consci_religios <- ifelse(gss$consci_cat_binary=="Percent 'A great deal' of trust in scientists", "Religiosity", 0)
gss$consci_educ <- ifelse(gss$consci_cat_binary=="Percent 'A great deal' of trust in scientists", "Education", 0)
gss$consci_income <- ifelse(gss$consci_cat_binary=="Percent 'A great deal' of trust in scientists", "Income Distribution Percentile", 0)
gss$consci_religion <- ifelse(gss$consci_cat_binary=="Percent 'A great deal' of trust in scientists", "Religion", 0)
gss$consci_parents_educ <- ifelse(gss$consci_cat_binary=="Percent 'A great deal' of trust in scientists", "Parents' Education", 0)
gss$consci_class <- ifelse(gss$consci_cat_binary=="Percent 'A great deal' of trust in scientists", "Class", 0)
gss$consci_norelig <- ifelse(gss$consci_cat_binary=="Percent 'A great deal' of trust in scientists", "Religion: None", 0)


p1a <- crosstab_3way(df = gss, x = sex, y = consci_gender, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, linetype = sex, col = sex)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(consci_gender)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p2a <- crosstab_3way(df = gss, x = race_binary, y = consci_race, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, linetype = race_binary, col = race_binary)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(consci_race)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p3a <- crosstab_3way(df = gss, x = rur_urb_binary, y = consci_urb, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, linetype = rur_urb_binary, col = rur_urb_binary)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(consci_urb)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p4a <- crosstab_3way(df = gss, x = religiosity, y = consci_religios, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, linetype = religiosity, col = religiosity)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(consci_religios)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p5a <- crosstab_3way(df = gss, x = education, y = consci_educ, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, linetype = education, col = education)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(consci_educ)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +  geom_point() +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p6a <- crosstab_3way(df = gss, x = class2, y = consci_class, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, linetype = class2, col = class2)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(consci_class)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") + 
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

Fig2 <- ggarrange(p1a, p2a, p3a, p4a, p5a, p6a, nrow = 2, ncol = 3, align = "v")
Fig2

# Figure 3
gss$democrat_cat <- ifelse(gss$democrat==1, "Percent Democrat", "Not Democrat")

gss$dems_gender <- ifelse(gss$democrat_cat=="Percent Democrat", "Gender", 0)
gss$dems_race <- ifelse(gss$democrat_cat=="Percent Democrat", "Race", 0)
gss$dems_urb <- ifelse(gss$democrat_cat=="Percent Democrat", "Rurality", 0)
gss$dems_religios <- ifelse(gss$democrat_cat=="Percent Democrat", "Religiosity", 0)
gss$dems_educ <- ifelse(gss$democrat_cat=="Percent Democrat", "Education", 0)
gss$dems_class <- ifelse(gss$democrat_cat=="Percent Democrat", "Class", 0)


p1b <- crosstab_3way(df = gss, x = sex, y = dems_gender, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = sex, linetype = sex)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(dems_gender)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p2b <- crosstab_3way(df = gss, x = race_binary, y = dems_race, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = race_binary, linetype = race_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(dems_race)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p3b <- crosstab_3way(df = gss, x = rur_urb_binary, y = dems_urb, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = rur_urb_binary, linetype = rur_urb_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(dems_urb)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p4b <- crosstab_3way(df = gss, x = religiosity, y = dems_religios, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = religiosity, linetype = religiosity)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(dems_religios)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p5b <- crosstab_3way(df = gss, x = education, y = dems_educ, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = education, linetype = education)) +
  scale_color_grey(start = 0.1, end = 0.5) +
  geom_smooth() +
  facet_wrap(facets = vars(dems_educ)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +  geom_point() +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p6b <- crosstab_3way(df = gss, x = class2, y = dems_class, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = class2, linetype = class2)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(dems_class)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") + 
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

Fig3 <- ggarrange(p1b, p2b, p3b, p4b, p5b, p6b, nrow = 2, ncol = 3, align = "v")
Fig3

# Figure 4
Fig4 <- crosstab_3way(df = gss, x = pid, y = consci_cat, z = year, weight = wtssps, format = "long") %>%
  ggplot(aes(year, pct/100, col = pid, linetype = pid)) +
  geom_smooth() + 
  scale_color_grey(start = 0.05, end = 0.6) +
  facet_wrap(facets = vars(consci_cat)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + ggtitle("") +
  theme(legend.position = "bottom")
Fig4

# Figure 5
#scientists
demsci <- crosstab_3way(df = gss, x = democrat, y = consci_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
demsci <- demsci[, c("year", "pct")]
colnames(demsci)[2] <- "dempct"

repsci <- crosstab_3way(df = gss, x = republican, y = consci_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
repsci <- repsci[, c("year", "pct")]
colnames(repsci)[2] <- "reppct"

scipol <- merge(demsci, repsci, by = "year")
scipol$gap <- scipol$dempct - scipol$reppct

p_scipol <- ggplot(data = scipol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Scientific community") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# military
demarmy <- crosstab_3way(df = gss, x = democrat, y = conarmy_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
demarmy <- demarmy[, c("year", "pct")]
colnames(demarmy)[2] <- "dempct"

reparmy <- crosstab_3way(df = gss, x = republican, y = conarmy_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
reparmy <- reparmy[, c("year", "pct")]
colnames(reparmy)[2] <- "reppct"

armypol <- merge(demarmy, reparmy, by = "year")
armypol$gap <- armypol$dempct - armypol$reppct

p_armypol <- ggplot(data = armypol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Military") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# major companies
dembus <- crosstab_3way(df = gss, x = democrat, y = conbus_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
dembus <- dembus[, c("year", "pct")]
colnames(dembus)[2] <- "dempct"

repbus <- crosstab_3way(df = gss, x = republican, y = conbus_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
repbus <- repbus[, c("year", "pct")]
colnames(repbus)[2] <- "reppct"

buspol <- merge(dembus, repbus, by = "year")
buspol$gap <- buspol$dempct - buspol$reppct

p_buspol <- ggplot(data = buspol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Major companies") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# organized religion
demclerg <- crosstab_3way(df = gss, x = democrat, y = conclerg_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
demclerg <- demclerg[, c("year", "pct")]
colnames(demclerg)[2] <- "dempct"

repclerg <- crosstab_3way(df = gss, x = republican, y = conclerg_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
repclerg <- repclerg[, c("year", "pct")]
colnames(repclerg)[2] <- "reppct"

clergpol <- merge(demclerg, repclerg, by = "year")
clergpol$gap <- clergpol$dempct - clergpol$reppct

p_clergpol <- ggplot(data = clergpol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Organized religion") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# education
demeduc <- crosstab_3way(df = gss, x = democrat, y = coneduc_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
demeduc <- demeduc[, c("year", "pct")]
colnames(demeduc)[2] <- "dempct"

repeduc <- crosstab_3way(df = gss, x = republican, y = coneduc_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
repeduc <- repeduc[, c("year", "pct")]
colnames(repeduc)[2] <- "reppct"

educpol <- merge(demeduc, repeduc, by = "year")
educpol$gap <- educpol$dempct - educpol$reppct

p_educpol <- ggplot(data = educpol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Education") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# banks and financial institutions
demfinan <- crosstab_3way(df = gss, x = democrat, y = confinan_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
demfinan <- demfinan[, c("year", "pct")]
colnames(demfinan)[2] <- "dempct"

repfinan <- crosstab_3way(df = gss, x = republican, y = confinan_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
repfinan <- repfinan[, c("year", "pct")]
colnames(repfinan)[2] <- "reppct"

finanpol <- merge(demfinan, repfinan, by = "year")
finanpol$gap <- finanpol$dempct - finanpol$reppct

p_finanpol <- ggplot(data = finanpol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Banks and financial institutions") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# supreme court
demjudge <- crosstab_3way(df = gss, x = democrat, y = conjudge_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
demjudge <- demjudge[, c("year", "pct")]
colnames(demjudge)[2] <- "dempct"

repjudge <- crosstab_3way(df = gss, x = republican, y = conjudge_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
repjudge <- repjudge[, c("year", "pct")]
colnames(repjudge)[2] <- "reppct"

judgepol <- merge(demjudge, repjudge, by = "year")
judgepol$gap <- judgepol$dempct - judgepol$reppct

p_judgepol <- ggplot(data = judgepol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Supreme Court") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# organized labor
demlabor <- crosstab_3way(df = gss, x = democrat, y = conlabor_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
demlabor <- demlabor[, c("year", "pct")]
colnames(demlabor)[2] <- "dempct"

replabor <- crosstab_3way(df = gss, x = republican, y = conlabor_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
replabor <- replabor[, c("year", "pct")]
colnames(replabor)[2] <- "reppct"

laborpol <- merge(demlabor, replabor, by = "year")
laborpol$gap <- laborpol$dempct - laborpol$reppct

p_laborpol <- ggplot(data = laborpol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Organized labor") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# Congress
demlegis <- crosstab_3way(df = gss, x = democrat, y = conlegis_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
demlegis <- demlegis[, c("year", "pct")]
colnames(demlegis)[2] <- "dempct"

replegis <- crosstab_3way(df = gss, x = republican, y = conlegis_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
replegis <- replegis[, c("year", "pct")]
colnames(replegis)[2] <- "reppct"

legispol <- merge(demlegis, replegis, by = "year")
legispol$gap <- legispol$dempct - legispol$reppct

p_legispol <- ggplot(data = legispol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Congress") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# Medicine
demmedic <- crosstab_3way(df = gss, x = democrat, y = conmedic_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
demmedic <- demmedic[, c("year", "pct")]
colnames(demmedic)[2] <- "dempct"

repmedic <- crosstab_3way(df = gss, x = republican, y = conmedic_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
repmedic <- repmedic[, c("year", "pct")]
colnames(repmedic)[2] <- "reppct"

medicpol <- merge(demmedic, repmedic, by = "year")
medicpol$gap <- medicpol$dempct - medicpol$reppct

p_medicpol <- ggplot(data = medicpol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Medicine") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# press
dempress <- crosstab_3way(df = gss, x = democrat, y = conpress_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
dempress <- dempress[, c("year", "pct")]
colnames(dempress)[2] <- "dempct"

reppress <- crosstab_3way(df = gss, x = republican, y = conpress_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
reppress <- reppress[, c("year", "pct")]
colnames(reppress)[2] <- "reppct"

presspol <- merge(dempress, reppress, by = "year")
presspol$gap <- presspol$dempct - presspol$reppct

p_presspol <- ggplot(data = presspol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Press") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# television
demtv <- crosstab_3way(df = gss, x = democrat, y = contv_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
demtv <- demtv[, c("year", "pct")]
colnames(demtv)[2] <- "dempct"

reptv <- crosstab_3way(df = gss, x = republican, y = contv_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
reptv <- reptv[, c("year", "pct")]
colnames(reptv)[2] <- "reppct"

tvpol <- merge(demtv, reptv, by = "year")
tvpol$gap <- tvpol$dempct - tvpol$reppct

p_tvpol <- ggplot(data = tvpol, aes(x = year, y = gap)) + 
  geom_smooth(color = "black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Television") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

Fig5 <- ggarrange(p_scipol, p_educpol, p_presspol, p_tvpol,
                      p_judgepol, p_medicpol, p_armypol, p_clergpol,
                      p_buspol, p_laborpol, p_finanpol, p_legispol,
                      nrow = 3, ncol = 4, align = "v")
Fig5 <- annotate_figure(Fig5,
                            left = text_grob("Democrats with Great Deal of Confidence -\nRepublicans with Great Deal of Confidence", size = 14, rot = 90))
Fig5

# Figure 6
fig6data <- read_dta("fig6data.dta")

scitrustmeans <- c(mean(fig6data$confsci[fig6data$rep==1]), mean(fig6data$confsci[fig6data$dem==1]),
                   mean(fig6data$pconfsci[fig6data$rep==1]), mean(fig6data$pconfsci[fig6data$dem==1]),
                   mean(fig6data$scidiff[fig6data$rep==1]), mean(fig6data$scidiff[fig6data$dem==1]))
parties <- c("Republicans", "Democrats", "Republicans", "Democrats", "Republicans", "Democrats")
cat <- c("Trust in Scientists", "Trust in Scientists", "Perception of Out-Party Trust in Scientists", "Perception of Out-Party Trust in Scientists", "Trust Differential", "Trust Differential")
sds <- c(sd(fig6data$confsci[fig6data$rep==1]), sd(fig6data$confsci[fig6data$dem==1]),
         sd(fig6data$pconfsci[fig6data$rep==1]), sd(fig6data$pconfsci[fig6data$dem==1]),
         sd(fig6data$scidiff[fig6data$rep==1]), sd(fig6data$scidiff[fig6data$dem==1]))
n <- c(742, 759, 742, 759, 742, 759)

df <- data.frame(scitrustmeans,parties,cat, sds, n)
df$se <- qnorm(0.975) * (df$sds/sqrt(df$n))
df$Lower_CI <- (df$scitrustmeans - df$se)
df$Upper_CI <- (df$scitrustmeans + df$se)

df$cat <- factor(df$cat, levels = c("Trust Differential", "Perception of Out-Party Trust in Scientists", "Trust in Scientists"))

Fig6 <- ggplot(df, aes(parties, scitrustmeans)) +
  geom_col(aes(fill = forcats::fct_rev(cat)), position = position_dodge(0.7), 
           color = "black", alpha = 0.7, width = 0.6) +
  geom_errorbar(aes(ymin = Lower_CI, ymax = Upper_CI, group = forcats::fct_rev(cat)), 
                position = position_dodge(0.7), width = 0.2) +
  scale_fill_manual(values = c("black", "gray30", "gray70")) +
  theme_minimal(base_size = 20) +
  theme(panel.grid.major.x = element_blank()) +
  labs(x = NULL, y = "Average Score (0 to 100 scale)") + ylim(0,100) +
  theme(legend.position="bottom") +
  theme(legend.title=element_blank()) +
  annotate("text", x = 0.765, y = 86.34, label = "76.34", color="black", size=6) +
  annotate("text", x = 1, y = 45.84, label = "35.84", color="black", size=6) +
  annotate("text", x = 1.235, y = 53.34, label = "43.34", color="black", size=6) +
  annotate("text", x = 1.765, y = 65.29, label = "55.29", color="black", size=6) +
  annotate("text", x = 2, y = 65.74, label = "55.74", color="black", size=6) +
  annotate("text", x = 2.235, y = 36.64, label = "26.64", color="black", size=6)

Fig6


# Figure 7. Regression Coefficients with Trust in Scientists and Researchers.
fig7data <- read_stata("fig7data")
fig7data_reps <- subset(fig7data, fig7data$rep==1)
fig7data_dems <- subset(fig7data, fig7data$dem==1)

fig7data$pid=case_when(
  fig7data$rep==1~"Republican",
  fig7data$dem==1~"Democrat",
)

fig7_ovr <- lm(scienctr ~ medicaltr + applied + basic + socialsci + climatetr + covidtr, data = fig7data)
fig7_dems <- lm(scienctr ~ medicaltr + applied + basic + socialsci + climatetr + covidtr, data = fig7data_dems)
fig7_reps <- lm(scienctr ~ medicaltr + applied + basic + socialsci + climatetr + covidtr, data = fig7data_reps)

Fig7 <- plot_summs(fig7_ovr,fig7_dems,fig7_reps,
           model.names=c("Overall","Democrats","Republicans"),
           colors = c("black", "gray40", "gray70"),
           coefs = c("Medical Scientists"="medicaltr",
                     "Applied Scientists"="applied",
                     "Natural Scientists"="basic",
                     "Social Scientists"="socialsci",
                     "Climate Scientists"="climatetr",
                     "COVID-19 Biochemists"="covidtr",
                     point.size = 8,
                     legend.title="", facet.label.pos="right")) +
  ggtitle("") + xlim(-0.2,0.6) + xlab("Coefficients") + ylab("Trust Target") +
  theme(legend.title=element_blank()) +
  theme(text = element_text(size = 16)) +
  theme(axis.text.x = element_text(size = 16), axis.text.y = element_text(size = 16))

Fig7


# Table 2
tab2data <- read_dta("tab2data")

tab2data$naturalim <- rowMeans(subset(tab2data, select=c(physicistim, chemim, biologim)), na.rm = TRUE)

mean(tab2data$scimethodim)
sd(tab2data$scimethodim)
mean(tab2data$scimethodim[tab2data$pidrep<3.5])
sd(tab2data$scimethodim[tab2data$pidrep<3.5])
mean(tab2data$scimethodim[tab2data$pidrep>4.5])
sd(tab2data$scimethodim[tab2data$pidrep>4.5])


mean(tab2data$doctorim)
sd(tab2data$doctorim)
mean(tab2data$doctorim[tab2data$pidrep<3.5])
sd(tab2data$doctorim[tab2data$pidrep<3.5])
mean(tab2data$doctorim[tab2data$pidrep>4.5])
sd(tab2data$doctorim[tab2data$pidrep>4.5])

t.test(tab2data$scimethodim, tab2data$doctorim, alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep<3.5], tab2data$doctorim[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep>4.5], tab2data$doctorim[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$engineerim)
sd(tab2data$engineerim)
mean(tab2data$engineerim[tab2data$pidrep<3.5])
sd(tab2data$engineerim[tab2data$pidrep<3.5])
mean(tab2data$engineerim[tab2data$pidrep>4.5])
sd(tab2data$engineerim[tab2data$pidrep>4.5])

t.test(tab2data$scimethodim, tab2data$engineerim, alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep<3.5], tab2data$engineerim[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep>4.5], tab2data$engineerim[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$naturalim)
sd(tab2data$naturalim)
mean(tab2data$naturalim[tab2data$pidrep<3.5])
sd(tab2data$naturalim[tab2data$pidrep<3.5])
mean(tab2data$naturalim[tab2data$pidrep>4.5])
sd(tab2data$naturalim[tab2data$pidrep>4.5])

t.test(tab2data$scimethodim, tab2data$naturalim, alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep<3.5], tab2data$naturalim[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep>4.5], tab2data$naturalim[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$climateim)
sd(tab2data$climateim)
mean(tab2data$climateim[tab2data$pidrep<3.5])
sd(tab2data$climateim[tab2data$pidrep<3.5])
mean(tab2data$climateim[tab2data$pidrep>4.5])
sd(tab2data$climateim[tab2data$pidrep>4.5])

t.test(tab2data$scimethodim, tab2data$climateim, alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep<3.5], tab2data$climateim[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep>4.5], tab2data$climateim[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$covidim)
sd(tab2data$covidim)
mean(tab2data$covidim[tab2data$pidrep<3.5])
sd(tab2data$covidim[tab2data$pidrep<3.5])
mean(tab2data$covidim[tab2data$pidrep>4.5])
sd(tab2data$covidim[tab2data$pidrep>4.5])

t.test(tab2data$scimethodim, tab2data$covidim, alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep<3.5], tab2data$covidim[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep>4.5], tab2data$covidim[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$professorim)
sd(tab2data$professorim)
mean(tab2data$professorim[tab2data$pidrep<3.5])
sd(tab2data$professorim[tab2data$pidrep<3.5])
mean(tab2data$professorim[tab2data$pidrep>4.5])
sd(tab2data$professorim[tab2data$pidrep>4.5])

t.test(tab2data$scimethodim, tab2data$professorim, alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep<3.5], tab2data$professorim[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep>4.5], tab2data$professorim[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$highschoolim)
sd(tab2data$highschoolim)
mean(tab2data$highschoolim[tab2data$pidrep<3.5])
sd(tab2data$highschoolim[tab2data$pidrep<3.5])
mean(tab2data$highschoolim[tab2data$pidrep>4.5])
sd(tab2data$highschoolim[tab2data$pidrep>4.5])

t.test(tab2data$scimethodim, tab2data$highschoolim, alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep<3.5], tab2data$highschoolim[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep>4.5], tab2data$highschoolim[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$popularim)
sd(tab2data$popularim)
mean(tab2data$popularim[tab2data$pidrep<3.5])
sd(tab2data$popularim[tab2data$pidrep<3.5])
mean(tab2data$popularim[tab2data$pidrep>4.5])
sd(tab2data$popularim[tab2data$pidrep>4.5])

t.test(tab2data$scimethodim, tab2data$popularim, alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep<3.5], tab2data$popularim[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$scimethodim[tab2data$pidrep>4.5], tab2data$popularim[tab2data$pidrep>4.5], alternative = "two.sided")

################
### Appendix ###
################

# Looking for Mode Effects
print(crosstab_3way(df = gss, x = mode, y = consci_binary, z = year, weight = wtssps, format = "long", remove = 0), n=98)

# Figure 1
# Predictors of confidence in scientists by decade
gss1970s <- subset(gss, gss$year<1979.5)
gss1980s <- subset(gss, gss$year>1979.5 & gss$year<1989.5)
gss1990s <- subset(gss, gss$year>1989.5 & gss$year<1999.5)
gss2000s <- subset(gss, gss$year>1999.5 & gss$year<2009.5)
gss2010s <- subset(gss, gss$year>2009.5 & gss$year<2019.5)

sci_70s <- lm(consci_binary ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss1970s)
sci_80s <- lm(consci_binary ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss1980s)
sci_90s <- lm(consci_binary ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss1990s)
sci_00s <- lm(consci_binary ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss2000s)
sci_10s <- lm(consci_binary ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss2010s)

plot_summs(sci_70s,sci_80s,sci_90s,sci_00s,sci_10s,
           model.names=c("1970s","1980s","1990s","2000s","2010s"),
           colors = c("gray0", "gray20", "gray40", "gray60", "gray80"),
           coefs = c("Gender: Male"="sexMale",
                     "Race: White"="race_binaryWhite",
                     "Area: Non-Rural"="rur_urb_binaryNon-Rural",
                     "Religion: Not Religious"="religiosityNot Religious",
                     "Education: College Degree"="educationCollege Degree",
                     "Class: Middle/Upper"="class2Middle/Upper",
                     point.size = 8,
                     legend.title="", facet.label.pos="right")) +
  ggtitle("") + xlim(-0.1, 0.25) + xlab("Estimate") +
  theme(legend.title=element_blank()) +
  theme(text = element_text(size = 15))

# Figure 2
# Looking at lowest confidence in science (7-20-24)
gss$consci_gender2 <- ifelse(gss$consci_cat_lowest=="Percent 'hardly any' of trust in scientists", "Gender", 0)
gss$consci_race2 <- ifelse(gss$consci_cat_lowest=="Percent 'hardly any' of trust in scientists", "Race", 0)
gss$consci_urb2 <- ifelse(gss$consci_cat_lowest=="Percent 'hardly any' of trust in scientists", "Rurality", 0)
gss$consci_religios2 <- ifelse(gss$consci_cat_lowest=="Percent 'hardly any' of trust in scientists", "Religiosity", 0)
gss$consci_educ2 <- ifelse(gss$consci_cat_lowest=="Percent 'hardly any' of trust in scientists", "Education", 0)
gss$consci_income2 <- ifelse(gss$consci_cat_lowest=="Percent 'hardly any' of trust in scientists", "Income Distribution Percentile", 0)
gss$consci_religion2 <- ifelse(gss$consci_cat_lowest=="Percent 'hardly any' of trust in scientists", "Religion", 0)
gss$consci_parents_educ2 <- ifelse(gss$consci_cat_lowest=="Percent 'hardly any' of trust in scientists", "Parents' Education", 0)
gss$consci_class2 <- ifelse(gss$consci_cat_lowest=="Percent 'hardly any' of trust in scientists", "Class", 0)
gss$consci_norelig2 <- ifelse(gss$consci_cat_lowest=="Percent 'hardly any' of trust in scientists", "Religion: None", 0)


p1b <- crosstab_3way(df = gss, x = sex, y = consci_gender2, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = sex, linetype = sex)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(consci_gender2)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,0.25)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p2b <- crosstab_3way(df = gss, x = race_binary, y = consci_race2, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = race_binary, linetype = race_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(consci_race2)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,0.25)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p3b <- crosstab_3way(df = gss, x = rur_urb_binary, y = consci_urb2, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = rur_urb_binary, linetype = rur_urb_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(consci_urb2)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,0.25)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p4b <- crosstab_3way(df = gss, x = religiosity, y = consci_religios2, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = religiosity, linetype = religiosity)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(consci_religios2)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,0.25)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p5b <- crosstab_3way(df = gss, x = education, y = consci_educ2, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = education, linetype = education)) +
  scale_color_grey(start = 0.1, end = 0.5) +
  geom_smooth() +
  facet_wrap(facets = vars(consci_educ2)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,0.25)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +  geom_point() +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p6b <- crosstab_3way(df = gss, x = class2, y = consci_class2, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = class2, linetype = class2)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(consci_class2)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,0.25)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") + 
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

fig_dem_lowest <- ggarrange(p1b, p2b, p3b, p4b, p5b, p6b, nrow = 2, ncol = 3, align = "v")
fig_dem_lowest


# Percent Identifying as Republicans by Party over Time
gss$republican_cat <- ifelse(gss$republican==1, "Percent republican", "Not republican")

gss$reps_gender <- ifelse(gss$republican_cat=="Percent republican", "Gender", 0)
gss$reps_race <- ifelse(gss$republican_cat=="Percent republican", "Race", 0)
gss$reps_urb <- ifelse(gss$republican_cat=="Percent republican", "Rurality", 0)
gss$reps_religios <- ifelse(gss$republican_cat=="Percent republican", "Religiosity", 0)
gss$reps_educ <- ifelse(gss$republican_cat=="Percent republican", "Education", 0)
gss$reps_class <- ifelse(gss$republican_cat=="Percent republican", "Class", 0)


p1c <- crosstab_3way(df = gss, x = sex, y = reps_gender, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = sex, linetype = sex)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(reps_gender)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p2c <- crosstab_3way(df = gss, x = race_binary, y = reps_race, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = race_binary, linetype = race_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(reps_race)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p3c <- crosstab_3way(df = gss, x = rur_urb_binary, y = reps_urb, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = rur_urb_binary, linetype = rur_urb_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(reps_urb)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p4c <- crosstab_3way(df = gss, x = religiosity, y = reps_religios, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = religiosity, linetype = religiosity)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(reps_religios)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p5c <- crosstab_3way(df = gss, x = education, y = reps_educ, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = education, linetype = education)) +
  scale_color_grey(start = 0.1, end = 0.5) +
  geom_smooth() +
  facet_wrap(facets = vars(reps_educ)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +  geom_point() +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p6c <- crosstab_3way(df = gss, x = class2, y = reps_class, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = class2, linetype = class2)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(reps_class)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") + 
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

fig_republican <- ggarrange(p1c, p2c, p3c, p4c, p5c, p6c, nrow = 2, ncol = 3, align = "v")
fig_republican <- annotate_figure(fig_republican,
                                  top=text_grob("Percent Republican by Group", size = 18))
fig_republican


# Percent Identifying as Independent over Time
gss$independent_cat <- ifelse(gss$independent==1, "Percent independent", "Not independent")

gss$inds_gender <- ifelse(gss$independent_cat=="Percent independent", "Gender", 0)
gss$inds_race <- ifelse(gss$independent_cat=="Percent independent", "Race", 0)
gss$inds_urb <- ifelse(gss$independent_cat=="Percent independent", "Rurality", 0)
gss$inds_religios <- ifelse(gss$independent_cat=="Percent independent", "Religiosity", 0)
gss$inds_educ <- ifelse(gss$independent_cat=="Percent independent", "Education", 0)
gss$inds_class <- ifelse(gss$independent_cat=="Percent independent", "Class", 0)


p1d <- crosstab_3way(df = gss, x = sex, y = inds_gender, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = sex, linetype = sex)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(inds_gender)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p2d <- crosstab_3way(df = gss, x = race_binary, y = inds_race, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = race_binary, linetype = race_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(inds_race)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p3d <- crosstab_3way(df = gss, x = rur_urb_binary, y = inds_urb, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = rur_urb_binary, linetype = rur_urb_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(inds_urb)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p4d <- crosstab_3way(df = gss, x = religiosity, y = inds_religios, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = religiosity, linetype = religiosity)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(inds_religios)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p5d <- crosstab_3way(df = gss, x = education, y = inds_educ, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = education, linetype = education)) +
  scale_color_grey(start = 0.1, end = 0.5) +
  geom_smooth() +
  facet_wrap(facets = vars(inds_educ)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +  geom_point() +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p6d <- crosstab_3way(df = gss, x = class2, y = inds_class, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = class2, linetype = class2)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(inds_class)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") + 
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

fig_independent <- ggarrange(p1d, p2d, p3d, p4d, p5d, p6d, nrow = 2, ncol = 3, align = "v")
fig_independent <- annotate_figure(fig_independent,
                                  top=text_grob("Percent Independent by Group", size = 18))
fig_independent

# Predictors of being a Democrat by decade
dem_70s <- lm(democrat ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss1970s)
dem_80s <- lm(democrat ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss1980s)
dem_90s <- lm(democrat ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss1990s)
dem_00s <- lm(democrat ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss2000s)
dem_10s <- lm(democrat ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss2010s)

plot_summs(dem_70s,dem_80s,dem_90s,dem_00s,dem_10s,
           model.names=c("1970s","1980s","1990s","2000s","2010s"),
           colors = c("gray0", "gray20", "gray40", "gray60", "gray80"),
           coefs = c("Gender: Male"="sexMale",
                     "Race: White"="race_binaryWhite",
                     "Area: Non-Rural"="rur_urb_binaryNon-Rural",
                     "Religion: Not Religious"="religiosityNot Religious",
                     "Education: College Degree"="educationCollege Degree",
                     "Class: Middle/Upper"="class2Middle/Upper",
                     point.size = 8,
                     legend.title="", facet.label.pos="right")) +
  ggtitle("") + xlim(-0.4, 0.2) + xlab("Estimate") +
  theme(legend.title=element_blank()) +
  theme(text = element_text(size = 15))


# Results With Independent Leaners Treated as Independents Rather than Partisans 
# Democrats
gss$democrat_noleaners_cat <- ifelse(gss$democrat_noleaners==1, "Percent Democrat", "Not Democrat")

gss$dems_noleaners_gender <- ifelse(gss$democrat_noleaners_cat=="Percent Democrat", "Gender", 0)
gss$dems_noleaners_race <- ifelse(gss$democrat_noleaners_cat=="Percent Democrat", "Race", 0)
gss$dems_noleaners_urb <- ifelse(gss$democrat_noleaners_cat=="Percent Democrat", "Rurality", 0)
gss$dems_noleaners_religios <- ifelse(gss$democrat_noleaners_cat=="Percent Democrat", "Religiosity", 0)
gss$dems_noleaners_educ <- ifelse(gss$democrat_noleaners_cat=="Percent Democrat", "Education", 0)
gss$dems_noleaners_class <- ifelse(gss$democrat_noleaners_cat=="Percent Democrat", "Class", 0)


p1a_app <- crosstab_3way(df = gss, x = sex, y = dems_noleaners_gender, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = sex, linetype = sex)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(dems_noleaners_gender)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p2a_app <- crosstab_3way(df = gss, x = race_binary, y = dems_noleaners_race, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = race_binary, linetype = race_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(dems_noleaners_race)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p3a_app <- crosstab_3way(df = gss, x = rur_urb_binary, y = dems_noleaners_urb, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = rur_urb_binary, linetype = rur_urb_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(dems_noleaners_urb)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p4a_app <- crosstab_3way(df = gss, x = religiosity, y = dems_noleaners_religios, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = religiosity, linetype = religiosity)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(dems_noleaners_religios)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p5a_app <- crosstab_3way(df = gss, x = education, y = dems_noleaners_educ, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = education, linetype = education)) +
  scale_color_grey(start = 0.1, end = 0.5) +
  geom_smooth() +
  facet_wrap(facets = vars(dems_noleaners_educ)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +  geom_point() +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p6a_app <- crosstab_3way(df = gss, x = class2, y = dems_noleaners_class, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = class2, linetype = class2)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(dems_noleaners_class)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") + 
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

ggarrange(p1a_app, p2a_app, p3a_app, p4a_app, p5a_app, p6a_app, nrow = 2, ncol = 3, align = "v")


# Republicans
gss$republican_noleaners_cat <- ifelse(gss$republican_noleaners==1, "Percent Republican", "Not Republican")

gss$reps_noleaners_gender <- ifelse(gss$republican_noleaners_cat=="Percent Republican", "Gender", 0)
gss$reps_noleaners_race <- ifelse(gss$republican_noleaners_cat=="Percent Republican", "Race", 0)
gss$reps_noleaners_urb <- ifelse(gss$republican_noleaners_cat=="Percent Republican", "Rurality", 0)
gss$reps_noleaners_religios <- ifelse(gss$republican_noleaners_cat=="Percent Republican", "Religiosity", 0)
gss$reps_noleaners_educ <- ifelse(gss$republican_noleaners_cat=="Percent Republican", "Education", 0)
gss$reps_noleaners_class <- ifelse(gss$republican_noleaners_cat=="Percent Republican", "Class", 0)


p1b_app <- crosstab_3way(df = gss, x = sex, y = reps_noleaners_gender, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = sex, linetype = sex)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(reps_noleaners_gender)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p2b_app <- crosstab_3way(df = gss, x = race_binary, y = reps_noleaners_race, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = race_binary, linetype = race_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(reps_noleaners_race)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p3b_app <- crosstab_3way(df = gss, x = rur_urb_binary, y = reps_noleaners_urb, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = rur_urb_binary, linetype = rur_urb_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(reps_noleaners_urb)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p4b_app <- crosstab_3way(df = gss, x = religiosity, y = reps_noleaners_religios, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = religiosity, linetype = religiosity)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(reps_noleaners_religios)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p5b_app <- crosstab_3way(df = gss, x = education, y = reps_noleaners_educ, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = education, linetype = education)) +
  scale_color_grey(start = 0.1, end = 0.5) +
  geom_smooth() +
  facet_wrap(facets = vars(reps_noleaners_educ)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +  geom_point() +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p6b_app <- crosstab_3way(df = gss, x = class2, y = reps_noleaners_class, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = class2, linetype = class2)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(reps_noleaners_class)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") + 
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

ggarrange(p1b_app, p2b_app, p3b_app, p4b_app, p5b_app, p6b_app, nrow = 2, ncol = 3, align = "v")


# independents
gss$independent_noleaners_cat <- ifelse(gss$independent_noleaners==1, "Percent Independent", "Not independent")

gss$inds_noleaners_gender <- ifelse(gss$independent_noleaners_cat=="Percent Independent", "Gender", 0)
gss$inds_noleaners_race <- ifelse(gss$independent_noleaners_cat=="Percent Independent", "Race", 0)
gss$inds_noleaners_urb <- ifelse(gss$independent_noleaners_cat=="Percent Independent", "Rurality", 0)
gss$inds_noleaners_religios <- ifelse(gss$independent_noleaners_cat=="Percent Independent", "Religiosity", 0)
gss$inds_noleaners_educ <- ifelse(gss$independent_noleaners_cat=="Percent Independent", "Education", 0)
gss$inds_noleaners_class <- ifelse(gss$independent_noleaners_cat=="Percent Independent", "Class", 0)


p1c_app <- crosstab_3way(df = gss, x = sex, y = inds_noleaners_gender, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = sex, linetype = sex)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(inds_noleaners_gender)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p2c_app <- crosstab_3way(df = gss, x = race_binary, y = inds_noleaners_race, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = race_binary, linetype = race_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(inds_noleaners_race)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p3c_app <- crosstab_3way(df = gss, x = rur_urb_binary, y = inds_noleaners_urb, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = rur_urb_binary, linetype = rur_urb_binary)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(inds_noleaners_urb)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p4c_app <- crosstab_3way(df = gss, x = religiosity, y = inds_noleaners_religios, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = religiosity, linetype = religiosity)) +
  geom_smooth() +  geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(inds_noleaners_religios)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p5c_app <- crosstab_3way(df = gss, x = education, y = inds_noleaners_educ, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = education, linetype = education)) +
  scale_color_grey(start = 0.1, end = 0.5) +
  geom_smooth() +
  facet_wrap(facets = vars(inds_noleaners_educ)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") +  geom_point() +
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p6c_app <- crosstab_3way(df = gss, x = class2, y = inds_noleaners_class, z = year, weight = wtssps, format = "long", remove = c("0")) %>%
  ggplot(aes(year, pct/100, col = class2, linetype = class2)) +
  geom_smooth() + geom_point() +
  scale_color_grey(start = 0.1, end = 0.5) +
  facet_wrap(facets = vars(inds_noleaners_class)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + 
  ggtitle("") + 
  theme(legend.position = "bottom") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

ggarrange(p1c_app, p2c_app, p3c_app, p4c_app, p5c_app, p6c_app, nrow = 2, ncol = 3, align = "v")

# Regression of Democrat Affiliation (Without Leaners as Partisans) By Decade 
gss1970s <- subset(gss, gss$year<1979.5)
gss1980s <- subset(gss, gss$year>1979.5 & gss$year<1989.5)
gss1990s <- subset(gss, gss$year>1989.5 & gss$year<1999.5)
gss2000s <- subset(gss, gss$year>1999.5 & gss$year<2009.5)
gss2010s <- subset(gss, gss$year>2009.5 & gss$year<2018.5)

dem_70s_app <- lm(democrat_noleaners ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss1970s)
dem_80s_app <- lm(democrat_noleaners ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss1980s)
dem_90s_app <- lm(democrat_noleaners ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss1990s)
dem_00s_app <- lm(democrat_noleaners ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss2000s)
dem_10s_app <- lm(democrat_noleaners ~ sex + race_binary + rur_urb_binary + religiosity + education + class2, data = gss2010s)

plot_summs(dem_70s_app,dem_80s_app,dem_90s_app,dem_00s_app,dem_10s_app,
           model.names=c("1970s","1980s","1990s","2000s","2010s"),
           colors = c("gray0", "gray20", "gray40", "gray60", "gray80"),
           coefs = c("Gender: Male"="sexMale",
                     "Race: White"="race_binaryWhite",
                     "Area: Non-Rural"="rur_urb_binaryNon-Rural",
                     "Religion: Not Religious"="religiosityNot Religious",
                     "Education: College Degree"="educationCollege Degree",
                     "Class: Middle/Upper"="class2Middle/Upper",
                     point.size = 8,
                     legend.title="", facet.label.pos="right")) +
  ggtitle("") + xlim(-0.4, 0.2) + xlab("Estimate") +
  theme(legend.title=element_blank()) +
  theme(text = element_text(size = 15))


# Confidence in the Scientific Community By Party (Without Leaners as Partisans) Over Time
crosstab_3way(df = gss, x = partyid_without_leaners_as_partisans, y = consci_cat, z = year, weight = wtssps, format = "long") %>%
  ggplot(aes(year, pct/100, col = partyid_without_leaners_as_partisans, linetype = partyid_without_leaners_as_partisans)) + 
  geom_smooth() + 
  scale_color_grey(start = 0.05, end = 0.6) +
  facet_wrap(facets = vars(consci_cat)) + ylab("") +
  scale_y_continuous(labels = scales::percent_format(accuracy=1), limits = c(0,1)) + theme_bw(base_size = 14) + xlab("") +
  theme(legend.title=element_blank()) + ggtitle("") +
  theme(legend.position = "bottom")


# Confidence in Scientists by Party (with leaners as partisans)
df <- crosstab_3way(df = gss, x = pid, y = consci_binary, z = year, weight = wtssps, remove=0, format = "long")

mean(df$pct[df$pid=="Democrat" & df$year==1973])
mean(df$pct[df$pid=="Republican" & df$year==1973])

mean(df$pct[df$pid=="Democrat" & df$year==2000])
mean(df$pct[df$pid=="Republican" & df$year==2000])

mean(df$pct[df$pid=="Democrat" & df$year==2018])
mean(df$pct[df$pid=="Republican" & df$year==2018])

mean(df$pct[df$pid=="Democrat" & df$year==2021])
mean(df$pct[df$pid=="Republican" & df$year==2021])

mean(df$pct[df$pid=="Democrat" & df$year==2022])
mean(df$pct[df$pid=="Republican" & df$year==2022])

mean(df$pct[df$pid=="Democrat" & df$year==2024])
mean(df$pct[df$pid=="Republican" & df$year==2024])



# Confidence in Scientists by Party (without leaners as partisans)
df <- crosstab_3way(df = gss, x = partyid_without_leaners_as_partisans, y = consci_binary, z = year, weight = wtssps, remove=0, format = "long")

mean(df$pct[df$partyid_without_leaners_as_partisans=="Democrat" & df$year==1973])
mean(df$pct[df$partyid_without_leaners_as_partisans=="Republican" & df$year==1973])

mean(df$pct[df$partyid_without_leaners_as_partisans=="Democrat" & df$year==2000])
mean(df$pct[df$partyid_without_leaners_as_partisans=="Republican" & df$year==2000])

mean(df$pct[df$partyid_without_leaners_as_partisans=="Democrat" & df$year==2018])
mean(df$pct[df$partyid_without_leaners_as_partisans=="Republican" & df$year==2018])

mean(df$pct[df$partyid_without_leaners_as_partisans=="Democrat" & df$year==2021])
mean(df$pct[df$partyid_without_leaners_as_partisans=="Republican" & df$year==2021])

mean(df$pct[df$partyid_without_leaners_as_partisans=="Democrat" & df$year==2022])
mean(df$pct[df$partyid_without_leaners_as_partisans=="Republican" & df$year==2022])

mean(df$pct[df$partyid_without_leaners_as_partisans=="Democrat" & df$year==2024])
mean(df$pct[df$partyid_without_leaners_as_partisans=="Republican" & df$year==2024])



# Partisan Gaps (Without Leaners as Partisans) in Confidence in Institutions over Time
#scientists
demsci <- crosstab_3way(df = gss, x = democrat_noleaners, y = consci_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
demsci <- demsci[, c("year", "pct")]
colnames(demsci)[2] <- "dempct"

repsci <- crosstab_3way(df = gss, x = republican_noleaners, y = consci_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
repsci <- repsci[, c("year", "pct")]
colnames(repsci)[2] <- "reppct"

scipol <- merge(demsci, repsci, by = "year")
scipol$gap <- scipol$dempct - scipol$reppct

p_scipol <- ggplot(data = scipol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Scientific community") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# military
demarmy <- crosstab_3way(df = gss, x = democrat_noleaners, y = conarmy_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
demarmy <- demarmy[, c("year", "pct")]
colnames(demarmy)[2] <- "dempct"

reparmy <- crosstab_3way(df = gss, x = republican_noleaners, y = conarmy_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
reparmy <- reparmy[, c("year", "pct")]
colnames(reparmy)[2] <- "reppct"

armypol <- merge(demarmy, reparmy, by = "year")
armypol$gap <- armypol$dempct - armypol$reppct

p_armypol <- ggplot(data = armypol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Military") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# major companies
dembus <- crosstab_3way(df = gss, x = democrat_noleaners, y = conbus_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
dembus <- dembus[, c("year", "pct")]
colnames(dembus)[2] <- "dempct"

repbus <- crosstab_3way(df = gss, x = republican_noleaners, y = conbus_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
repbus <- repbus[, c("year", "pct")]
colnames(repbus)[2] <- "reppct"

buspol <- merge(dembus, repbus, by = "year")
buspol$gap <- buspol$dempct - buspol$reppct

p_buspol <- ggplot(data = buspol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Major companies") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# organized religion
demclerg <- crosstab_3way(df = gss, x = democrat_noleaners, y = conclerg_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
demclerg <- demclerg[, c("year", "pct")]
colnames(demclerg)[2] <- "dempct"

repclerg <- crosstab_3way(df = gss, x = republican_noleaners, y = conclerg_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
repclerg <- repclerg[, c("year", "pct")]
colnames(repclerg)[2] <- "reppct"

clergpol <- merge(demclerg, repclerg, by = "year")
clergpol$gap <- clergpol$dempct - clergpol$reppct

p_clergpol <- ggplot(data = clergpol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Organized religion") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# education
demeduc <- crosstab_3way(df = gss, x = democrat_noleaners, y = coneduc_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
demeduc <- demeduc[, c("year", "pct")]
colnames(demeduc)[2] <- "dempct"

repeduc <- crosstab_3way(df = gss, x = republican_noleaners, y = coneduc_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
repeduc <- repeduc[, c("year", "pct")]
colnames(repeduc)[2] <- "reppct"

educpol <- merge(demeduc, repeduc, by = "year")
educpol$gap <- educpol$dempct - educpol$reppct

p_educpol <- ggplot(data = educpol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Education") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# executive branch
demfed <- crosstab_3way(df = gss, x = democrat_noleaners, y = confed_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
demfed <- demfed[, c("year", "pct")]
colnames(demfed)[2] <- "dempct"

repfed <- crosstab_3way(df = gss, x = republican_noleaners, y = confed_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
repfed <- repfed[, c("year", "pct")]
colnames(repfed)[2] <- "reppct"

fedpol <- merge(demfed, repfed, by = "year")
fedpol$gap <- fedpol$dempct - fedpol$reppct

p_fedpol <- ggplot(data = fedpol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Partisan Gap in Executive Branch") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

p_fedpol <- annotate_figure(p_fedpol,
                            left = text_grob("democrat_noleanerss with Great Deal of Confidence -\nrepublican_noleanerss with Great Deal of Confidence", size = 14, rot = 90))

# banks and financial institutions
demfinan <- crosstab_3way(df = gss, x = democrat_noleaners, y = confinan_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
demfinan <- demfinan[, c("year", "pct")]
colnames(demfinan)[2] <- "dempct"

repfinan <- crosstab_3way(df = gss, x = republican_noleaners, y = confinan_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
repfinan <- repfinan[, c("year", "pct")]
colnames(repfinan)[2] <- "reppct"

finanpol <- merge(demfinan, repfinan, by = "year")
finanpol$gap <- finanpol$dempct - finanpol$reppct

p_finanpol <- ggplot(data = finanpol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Banks and financial institutions") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# supreme court
demjudge <- crosstab_3way(df = gss, x = democrat_noleaners, y = conjudge_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
demjudge <- demjudge[, c("year", "pct")]
colnames(demjudge)[2] <- "dempct"

repjudge <- crosstab_3way(df = gss, x = republican_noleaners, y = conjudge_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
repjudge <- repjudge[, c("year", "pct")]
colnames(repjudge)[2] <- "reppct"

judgepol <- merge(demjudge, repjudge, by = "year")
judgepol$gap <- judgepol$dempct - judgepol$reppct

p_judgepol <- ggplot(data = judgepol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Supreme Court") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# organized labor
demlabor <- crosstab_3way(df = gss, x = democrat_noleaners, y = conlabor_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
demlabor <- demlabor[, c("year", "pct")]
colnames(demlabor)[2] <- "dempct"

replabor <- crosstab_3way(df = gss, x = republican_noleaners, y = conlabor_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
replabor <- replabor[, c("year", "pct")]
colnames(replabor)[2] <- "reppct"

laborpol <- merge(demlabor, replabor, by = "year")
laborpol$gap <- laborpol$dempct - laborpol$reppct

p_laborpol <- ggplot(data = laborpol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Organized labor") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# Congress
demlegis <- crosstab_3way(df = gss, x = democrat_noleaners, y = conlegis_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
demlegis <- demlegis[, c("year", "pct")]
colnames(demlegis)[2] <- "dempct"

replegis <- crosstab_3way(df = gss, x = republican_noleaners, y = conlegis_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
replegis <- replegis[, c("year", "pct")]
colnames(replegis)[2] <- "reppct"

legispol <- merge(demlegis, replegis, by = "year")
legispol$gap <- legispol$dempct - legispol$reppct

p_legispol <- ggplot(data = legispol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Congress") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# Medicine
demmedic <- crosstab_3way(df = gss, x = democrat_noleaners, y = conmedic_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
demmedic <- demmedic[, c("year", "pct")]
colnames(demmedic)[2] <- "dempct"

repmedic <- crosstab_3way(df = gss, x = republican_noleaners, y = conmedic_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
repmedic <- repmedic[, c("year", "pct")]
colnames(repmedic)[2] <- "reppct"

medicpol <- merge(demmedic, repmedic, by = "year")
medicpol$gap <- medicpol$dempct - medicpol$reppct

p_medicpol <- ggplot(data = medicpol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Medicine") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# press
dempress <- crosstab_3way(df = gss, x = democrat_noleaners, y = conpress_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
dempress <- dempress[, c("year", "pct")]
colnames(dempress)[2] <- "dempct"

reppress <- crosstab_3way(df = gss, x = republican_noleaners, y = conpress_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
reppress <- reppress[, c("year", "pct")]
colnames(reppress)[2] <- "reppct"

presspol <- merge(dempress, reppress, by = "year")
presspol$gap <- presspol$dempct - presspol$reppct

p_presspol <- ggplot(data = presspol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Press") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

# television
demtv <- crosstab_3way(df = gss, x = democrat_noleaners, y = contv_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
demtv <- demtv[, c("year", "pct")]
colnames(demtv)[2] <- "dempct"

reptv <- crosstab_3way(df = gss, x = republican_noleaners, y = contv_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
reptv <- reptv[, c("year", "pct")]
colnames(reptv)[2] <- "reppct"

tvpol <- merge(demtv, reptv, by = "year")
tvpol$gap <- tvpol$dempct - tvpol$reppct

p_tvpol <- ggplot(data = tvpol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("Television") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

FigPartisanGap_app <- ggarrange(p_scipol, p_educpol, p_presspol, p_tvpol,
                  p_judgepol, p_medicpol, p_armypol, p_clergpol,
                  p_buspol, p_laborpol, p_finanpol, p_legispol,
                  nrow = 3, ncol = 4, align = "v")
FigPartisanGap_app <- annotate_figure(FigPartisanGap_app,
                        left = text_grob("Democrats with Great Deal of Confidence -\nRepublicans with Great Deal of Confidence", size = 14, rot = 90))
FigPartisanGap_app

# Executive branch plots
demfed <- crosstab_3way(df = gss, x = democrat, y = confed_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
demfed <- demfed[, c("year", "pct")]
colnames(demfed)[2] <- "dempct"

repfed <- crosstab_3way(df = gss, x = republican, y = confed_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
repfed <- repfed[, c("year", "pct")]
colnames(repfed)[2] <- "reppct"

fedpol <- merge(demfed, repfed, by = "year")
fedpol$gap <- fedpol$dempct - fedpol$reppct

p_fedpol <- ggplot(data = fedpol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

annotate_figure(p_fedpol,
                left = text_grob("Democrats with Great Deal of Confidence -\nRepublicans with Great Deal of Confidence", size = 14, rot = 90))

# Executive branch plot without leaners as partisans 
demfed <- crosstab_3way(df = gss, x = democrat_noleaners, y = confed_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
demfed <- demfed[, c("year", "pct")]
colnames(demfed)[2] <- "dempct"

repfed <- crosstab_3way(df = gss, x = republican_noleaners, y = confed_binary, z = year, weight = wtssps, format = "long", remove = c("0"))
repfed <- repfed[, c("year", "pct")]
colnames(repfed)[2] <- "reppct"

fedpol <- merge(demfed, repfed, by = "year")
fedpol$gap <- fedpol$dempct - fedpol$reppct

p_fedpol <- ggplot(data = fedpol, aes(x = year, y = gap)) + 
  geom_smooth(color="black") +
  geom_point() + xlab("") + ylab("") + ylim(-40,40) +
  theme_bw(base_size = 12) + ggtitle("") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=0, linetype="dashed")

annotate_figure(p_fedpol,
                left = text_grob("Democrats with Great Deal of Confidence -\nRepublicans with Great Deal of Confidence", size = 14, rot = 90))



# Perceptions of Institutions Data

# Trust Differentials table in appendix
# scientific community
mean(fig6data$scidiff)
sd(fig6data$scidiff)
mean(fig6data$scidiff[fig6data$dem==1])
sd(fig6data$scidiff[fig6data$dem==1])
mean(fig6data$scidiff[fig6data$rep==1])
sd(fig6data$scidiff[fig6data$rep==1])

# 
mean(fig6data$educdiff)
sd(fig6data$educdiff)
mean(fig6data$educdiff[fig6data$dem==1])
sd(fig6data$educdiff[fig6data$dem==1])
mean(fig6data$educdiff[fig6data$rep==1])
sd(fig6data$educdiff[fig6data$rep==1])

# 
mean(fig6data$mediadiff)
sd(fig6data$mediadiff)
mean(fig6data$mediadiff[fig6data$dem==1])
sd(fig6data$mediadiff[fig6data$dem==1])
mean(fig6data$mediadiff[fig6data$rep==1])
sd(fig6data$mediadiff[fig6data$rep==1])

# 
mean(fig6data$socmeddiff)
sd(fig6data$socmeddiff)
mean(fig6data$socmeddiff[fig6data$dem==1])
sd(fig6data$socmeddiff[fig6data$dem==1])
mean(fig6data$socmeddiff[fig6data$rep==1])
sd(fig6data$socmeddiff[fig6data$rep==1])

# 
mean(fig6data$militarydiff)
sd(fig6data$militarydiff)
mean(fig6data$militarydiff[fig6data$dem==1])
sd(fig6data$militarydiff[fig6data$dem==1])
mean(fig6data$militarydiff[fig6data$rep==1])
sd(fig6data$militarydiff[fig6data$rep==1])

# 
mean(fig6data$policediff)
sd(fig6data$policediff)
mean(fig6data$policediff[fig6data$dem==1])
sd(fig6data$policediff[fig6data$dem==1])
mean(fig6data$policediff[fig6data$rep==1])
sd(fig6data$policediff[fig6data$rep==1])

# 
mean(fig6data$docotorsdiff)
sd(fig6data$docotorsdiff)
mean(fig6data$docotorsdiff[fig6data$dem==1])
sd(fig6data$docotorsdiff[fig6data$dem==1])
mean(fig6data$docotorsdiff[fig6data$rep==1])
sd(fig6data$docotorsdiff[fig6data$rep==1])

# 
mean(fig6data$religdiff)
sd(fig6data$religdiff)
mean(fig6data$religdiff[fig6data$dem==1])
sd(fig6data$religdiff[fig6data$dem==1])
mean(fig6data$religdiff[fig6data$rep==1])
sd(fig6data$religdiff[fig6data$rep==1])

# 
mean(fig6data$labordiff)
sd(fig6data$labordiff)
mean(fig6data$labordiff[fig6data$dem==1])
sd(fig6data$labordiff[fig6data$dem==1])
mean(fig6data$labordiff[fig6data$rep==1])
sd(fig6data$labordiff[fig6data$rep==1])

# 
mean(fig6data$companiesdiff)
sd(fig6data$companiesdiff)
mean(fig6data$companiesdiff[fig6data$dem==1])
sd(fig6data$companiesdiff[fig6data$dem==1])
mean(fig6data$companiesdiff[fig6data$rep==1])
sd(fig6data$companiesdiff[fig6data$rep==1])

# 
mean(fig6data$courtdiff)
sd(fig6data$courtdiff)
mean(fig6data$courtdiff[fig6data$dem==1])
sd(fig6data$courtdiff[fig6data$dem==1])
mean(fig6data$courtdiff[fig6data$rep==1])
sd(fig6data$courtdiff[fig6data$rep==1])

# 
mean(fig6data$congressdiff)
sd(fig6data$congressdiff)
mean(fig6data$congressdiff[fig6data$dem==1])
sd(fig6data$congressdiff[fig6data$dem==1])
mean(fig6data$congressdiff[fig6data$rep==1])
sd(fig6data$congressdiff[fig6data$rep==1])

# 
mean(fig6data$presidentdiff)
sd(fig6data$presidentdiff)
mean(fig6data$presidentdiff[fig6data$dem==1])
sd(fig6data$presidentdiff[fig6data$dem==1])
mean(fig6data$presidentdiff[fig6data$rep==1])
sd(fig6data$presidentdiff[fig6data$rep==1])

# Who are Scientists Data

# Tables Underlying Coefficient Plots
# Figure 7
summary(fig7_ovr)
summary(fig7_dems)
summary(fig7_reps)

stargazer::stargazer(fig7_ovr, fig7_dems, fig7_reps,
                     column.labels=c("All Respondents", "Democrats", "Republicans"),
                     dep.var.labels = "General Trust in Scientists",
                     covariate.labels=c("Constant", "Medical Scientists", "Applied Scientists", "Natural Scientists",
                                        "Social Scientists", "Climate Scientists", "COVID-19 Biochemists"),
                     summary = TRUE, align=TRUE,
                     type="html",
                     font.size="scriptsize", 
                     table.placement="H",
                     column.sep.width = "1pt", 
                     dep.var.labels.include = T,
                     intercept.bottom = FALSE,
                     out = "C:/Users/jonsc/Downloads/test.html", report=('vcsp'),
                     digits = 3,  no.space=TRUE)


# Trust Scores
mean(tab2data$scienctr)
sd(tab2data$scienctr)
mean(tab2data$scienctr[tab2data$pidrep<3.5])
sd(tab2data$scienctr[tab2data$pidrep<3.5])
mean(tab2data$scienctr[tab2data$pidrep>4.5])
sd(tab2data$scienctr[tab2data$pidrep>4.5])

mean(tab2data$medicaltr)
sd(tab2data$medicaltr)
mean(tab2data$medicaltr[tab2data$pidrep<3.5])
sd(tab2data$medicaltr[tab2data$pidrep<3.5])
mean(tab2data$medicaltr[tab2data$pidrep>4.5])
sd(tab2data$medicaltr[tab2data$pidrep>4.5])

t.test(tab2data$medicaltr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$medicaltr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$medicaltr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

tab2data$appliedtr <- rowMeans(subset(tab2data, select=c(aerotr, civiltr, computertr)), na.rm = TRUE)

mean(tab2data$appliedtr)
sd(tab2data$appliedtr)
mean(tab2data$appliedtr[tab2data$pidrep<3.5])
sd(tab2data$appliedtr[tab2data$pidrep<3.5])
mean(tab2data$appliedtr[tab2data$pidrep>4.5])
sd(tab2data$appliedtr[tab2data$pidrep>4.5])

t.test(tab2data$appliedtr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$appliedtr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$appliedtr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$aerotr)
sd(tab2data$aerotr)
mean(tab2data$aerotr[tab2data$pidrep<3.5])
sd(tab2data$aerotr[tab2data$pidrep<3.5])
mean(tab2data$aerotr[tab2data$pidrep>4.5])
sd(tab2data$aerotr[tab2data$pidrep>4.5])

t.test(tab2data$aerotr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$aerotr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$aerotr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$civiltr)
sd(tab2data$civiltr)
mean(tab2data$civiltr[tab2data$pidrep<3.5])
sd(tab2data$civiltr[tab2data$pidrep<3.5])
mean(tab2data$civiltr[tab2data$pidrep>4.5])
sd(tab2data$civiltr[tab2data$pidrep>4.5])

t.test(tab2data$civiltr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$civiltr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$civiltr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$computertr)
sd(tab2data$computertr)
mean(tab2data$computertr[tab2data$pidrep<3.5])
sd(tab2data$computertr[tab2data$pidrep<3.5])
mean(tab2data$computertr[tab2data$pidrep>4.5])
sd(tab2data$computertr[tab2data$pidrep>4.5])

t.test(tab2data$computertr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$computertr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$computertr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

tab2data$naturaltr <- rowMeans(subset(tab2data, select=c(physictr, chemisttr, biologisttr)), na.rm = TRUE)

mean(tab2data$naturaltr)
sd(tab2data$naturaltr)
mean(tab2data$naturaltr[tab2data$pidrep<3.5])
sd(tab2data$naturaltr[tab2data$pidrep<3.5])
mean(tab2data$naturaltr[tab2data$pidrep>4.5])
sd(tab2data$naturaltr[tab2data$pidrep>4.5])

t.test(tab2data$naturaltr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$naturaltr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$naturaltr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$biologisttr)
sd(tab2data$biologisttr)
mean(tab2data$biologisttr[tab2data$pidrep<3.5])
sd(tab2data$biologisttr[tab2data$pidrep<3.5])
mean(tab2data$biologisttr[tab2data$pidrep>4.5])
sd(tab2data$biologisttr[tab2data$pidrep>4.5])

t.test(tab2data$biologisttr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$biologisttr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$biologisttr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$chemisttr)
sd(tab2data$chemisttr)
mean(tab2data$chemisttr[tab2data$pidrep<3.5])
sd(tab2data$chemisttr[tab2data$pidrep<3.5])
mean(tab2data$chemisttr[tab2data$pidrep>4.5])
sd(tab2data$chemisttr[tab2data$pidrep>4.5])

t.test(tab2data$chemisttr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$chemisttr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$chemisttr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$physictr)
sd(tab2data$physictr)
mean(tab2data$physictr[tab2data$pidrep<3.5])
sd(tab2data$physictr[tab2data$pidrep<3.5])
mean(tab2data$physictr[tab2data$pidrep>4.5])
sd(tab2data$physictr[tab2data$pidrep>4.5])

t.test(tab2data$physictr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$physictr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$physictr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

tab2data$socscientiststr <- rowMeans(subset(tab2data, select=c(psychtr, econtr, poliscitr)), na.rm = TRUE)

mean(tab2data$socscientiststr)
sd(tab2data$socscientiststr)
mean(tab2data$socscientiststr[tab2data$pidrep<3.5])
sd(tab2data$socscientiststr[tab2data$pidrep<3.5])
mean(tab2data$socscientiststr[tab2data$pidrep>4.5])
sd(tab2data$socscientiststr[tab2data$pidrep>4.5])

t.test(tab2data$socscientiststr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$socscientiststr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$socscientiststr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$psychtr)
sd(tab2data$psychtr)
mean(tab2data$psychtr[tab2data$pidrep<3.5])
sd(tab2data$psychtr[tab2data$pidrep<3.5])
mean(tab2data$psychtr[tab2data$pidrep>4.5])
sd(tab2data$psychtr[tab2data$pidrep>4.5])

t.test(tab2data$psychtr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$psychtr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$psychtr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$econtr)
sd(tab2data$econtr)
mean(tab2data$econtr[tab2data$pidrep<3.5])
sd(tab2data$econtr[tab2data$pidrep<3.5])
mean(tab2data$econtr[tab2data$pidrep>4.5])
sd(tab2data$econtr[tab2data$pidrep>4.5])

t.test(tab2data$econtr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$econtr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$econtr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$poliscitr)
sd(tab2data$poliscitr)
mean(tab2data$poliscitr[tab2data$pidrep<3.5])
sd(tab2data$poliscitr[tab2data$pidrep<3.5])
mean(tab2data$poliscitr[tab2data$pidrep>4.5])
sd(tab2data$poliscitr[tab2data$pidrep>4.5])

t.test(tab2data$poliscitr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$poliscitr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$poliscitr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$climatetr)
sd(tab2data$climatetr)
mean(tab2data$climatetr[tab2data$pidrep<3.5])
sd(tab2data$climatetr[tab2data$pidrep<3.5])
mean(tab2data$climatetr[tab2data$pidrep>4.5])
sd(tab2data$climatetr[tab2data$pidrep>4.5])

t.test(tab2data$climatetr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$climatetr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$climatetr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$covidtr)
sd(tab2data$covidtr)
mean(tab2data$covidtr[tab2data$pidrep<3.5])
sd(tab2data$covidtr[tab2data$pidrep<3.5])
mean(tab2data$covidtr[tab2data$pidrep>4.5])
sd(tab2data$covidtr[tab2data$pidrep>4.5])

t.test(tab2data$covidtr, tab2data$scienctr, alternative = "two.sided")
t.test(tab2data$covidtr[tab2data$pidrep<3.5], tab2data$scienctr[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$covidtr[tab2data$pidrep>4.5], tab2data$scienctr[tab2data$pidrep>4.5], alternative = "two.sided")

# Perception Scores
mean(tab2data$scimethodim)
sd(tab2data$scimethodim)
mean(tab2data$scimethodim[tab2data$pidrep<3.5])
sd(tab2data$scimethodim[tab2data$pidrep<3.5])
mean(tab2data$scimethodim[tab2data$pidrep>4.5])
sd(tab2data$scimethodim[tab2data$pidrep>4.5])

mean(tab2data$biologim)
sd(tab2data$biologim)
mean(tab2data$biologim[tab2data$pidrep<3.5])
sd(tab2data$biologim[tab2data$pidrep<3.5])
mean(tab2data$biologim[tab2data$pidrep>4.5])
sd(tab2data$biologim[tab2data$pidrep>4.5])

t.test(tab2data$biologim, tab2data$scimethodim, alternative = "two.sided")
t.test(tab2data$biologim[tab2data$pidrep<3.5], tab2data$scimethodim[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$biologim[tab2data$pidrep>4.5], tab2data$scimethodim[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$chemim)
sd(tab2data$chemim)
mean(tab2data$chemim[tab2data$pidrep<3.5])
sd(tab2data$chemim[tab2data$pidrep<3.5])
mean(tab2data$chemim[tab2data$pidrep>4.5])
sd(tab2data$chemim[tab2data$pidrep>4.5])

t.test(tab2data$chemim, tab2data$scimethodim, alternative = "two.sided")
t.test(tab2data$chemim[tab2data$pidrep<3.5], tab2data$scimethodim[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$chemim[tab2data$pidrep>4.5], tab2data$scimethodim[tab2data$pidrep>4.5], alternative = "two.sided")

mean(tab2data$physicistim)
sd(tab2data$physicistim)
mean(tab2data$physicistim[tab2data$pidrep<3.5])
sd(tab2data$physicistim[tab2data$pidrep<3.5])
mean(tab2data$physicistim[tab2data$pidrep>4.5])
sd(tab2data$physicistim[tab2data$pidrep>4.5])

t.test(tab2data$physicistim, tab2data$scimethodim, alternative = "two.sided")
t.test(tab2data$physicistim[tab2data$pidrep<3.5], tab2data$scimethodim[tab2data$pidrep<3.5], alternative = "two.sided")
t.test(tab2data$physicistim[tab2data$pidrep>4.5], tab2data$scimethodim[tab2data$pidrep>4.5], alternative = "two.sided")
